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We present a simulation algorithm for Hamiltonian fermion lattice models. A guiding trial wave 
function is adaptively optimized during Monte Carlo evolution. We apply the method to the two 
dimensional Gross-Neveu model and analyze systematc errors in the study of ground state properties. 
We show that accurate measurements can be achieved by a proper extrapolation in the algorithm 
free parameters. 
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Lattice Field Theory is a constructive framework where 
non-perturbative properties of quantum models can be 
addressed both analytically and by numerical techniques. 
The main standing theoretical viewpoints are the tra- 
ditional Lagrangian approach jlj and the Hamiltonian 
formulation In the study of fermionic models, La- 
grangian simulations suffer the drawback of requiring 
Grassmann variables that are difficult to handle numer- 
ically and must be integrated out explicitly leading to 
large non-local determinants. Instead, in the Hamilto- 
nian approach, the treatment of Fermi anticommuting 
operators is straightforward. In particular, this holds in 
one spatial dimension where notoriously difficult sign- 
problems |3j are tame. 

Another important reason to resort to Hamiltonian 
methods is that they rely on powerful well founded Many- 
Body techniques Q . In particular, a direct analysis of the 
ground state structure is often feasible through a guid- 
ing trial wave function ||. This is an approximation 
to the exact ground state that can provide deep physi- 
cal insights about the model under consideration. Also, 
it plays a central role in the simulation algorithms and 
the quality of the results depends critically on its accu- 
racy |tJ . Usually, it contains a set of free parameters that 
deserve optimization by rather expensive variational cal- 
culations ||. 

Here, we present a Monte Carlo (MC) algorithm that 
includes automatic optimization of the trial wave func- 
tion by means of a non-linear feedback between state 
sampling and guiding. The MC core is based on a general 
stochastic representation of matrix evolution problems ||] 
and has been discussed in the specific case of the Hub- 
bard model jlO| . The adaptive optimization strategy has 
been already applied to Diffusion MC studies of purely 
bosonic models with continuous state space [||. 

In this Report, we focus on fermionic models and 
present an algorithm suitable for the study of Hamilto- 
nians acting on a finite-dimensional fully discrete state 
space. In fact, for a local fermion model discretized on 
a finite lattice, the Hamiltonian is a large sparse ma- 
trix H = {H ss i} SiS > e s, with S denoting the discrete 



state space. The ground state can be obtained by act- 
ing on a given initial state with the evolution semigroup 
O = {e~ tH } t >o in the t — > oo limit. For simplicity, we as- 
sume a non degenerate ground state, in the general case 
f2 projects onto the lowest eigenspace. 

To build a MC algorithm, we need a probabilistic rep- 
resentation of f2. For each pair s,s' G S such that 



s =/= s' and H s > s ^ we define T s 



We as- 



sume that all T s : s > (no sign-problem) and build a S- 
valued Markov stochastic process St by identifying r s / s 
as the rate for the transition s — > s'. Hence, the aver- 
age occupation P s (t) = E (5 s . St ) , with E (•) denoting the 
average with respect to St, obeys the Master Equation 

^(/?)=E sV s( r -' P s'- r ^)- 

Related to s t , we also define the real valued stochastic 
process W t = exp (- J Q * Lu St dtj , with uj s = J^s'es 
It can be shown that the weighted expectation value 
ip s (t) = E (<5 SjSt W t ) reconstructs fi: 



s'es 



with ip s (0) = Prob(so = s). Matrix elements of f2 can be 
identified with certain expectation values. In particular, 
the ground state energy Eq can be obtained by 

E(w st W t ) 



E 



lim 

t—>-\-oc 



(1) 



E (W t ) 

that gives Eq as the asymptotic average of u) s over real- 
izations of St with weight Wt, called walkers in the fol- 
lowing. The actual construction of the process is straigh- 
forward. A realization of St is a piece-wise constant map 
R — > S with isolated jumps at times t = to, t\, . . ., with 
to < t\ < t-z < ■ ■ ■. An algorithm to compute the triples 
{t n , s tn , W tn } is the following: 

1. We simply denote St n = s and define the set T s of 
target states connected to s: T s = {s',T s / s > 0}. 
We also define the total width T s = X) s 'eT ^s's- 

2. Extract r > with probability density p s (t) = 
r s e~ rsr . In other words, r = —J- log£ with £ uni- 
formly distributed in [0, 1]. 
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3. Extract a new state s' 6 T s with probability p s i = 

4. Define f n+ i = £„ + r, s tn+1 = s' and W tn+1 = 
W tn • e"^ T . 

The above algorithm is the explicit zero imaginary time 
limit of power algorithms OS . 

For a better performance, it is useful to introduce 
a trial state |$(a)) depending on some parameters a. 
The original Hamiltonian H is replaced by the isospec- 
tral H ss ,(a) = ^(^iJ^^^^a) with $ s (a) = (s|$(a)). 
The algorithm is unchanged (hermiticity of H (a) has not 
been assumed), but everything, in particular ui s , becomes 
a-dependent. In the ideal case when |$(a)) is the exact 
ground state, then ui s = E and the ground state energy 
is estimated by Eq. (Q) with zero fluctuations. 

As is well known, a naive implementation of Eq. (|l|) 
fails because the variance of the right hand side diverges 
as t -> +oo. A possible way out is Stochastic Recon- 
figuration (SR) An ensemble with a large fixed 
number K of walkers is introduced and a branching pro- 
cedure deletes walkers with low weight and makes copies 
of the ones with larger weight. In the end, we take the 
numerical limit K — » oo. If j3 is the time between two SR, 
then we denote the estimate of the ground state energy 
by Eq(P, K, a) where we do not write the dependence 
on physical parameters (lattice size, couplings). Usually, 
the dependence on a is quite strong and requires opti- 
mization to make |$(a)) the closest possible to the exact 
ground state. 

As we remarked, a possible way to optimize a is to 
minimize the fluctuations of uj St (a) Jl5[ |. To this aim, 
following the general ideas of p6| , we promote a to a se- 
quence {a n } and after each SR, we compute the variance 
of ui(a) over the K walkers, with their states kept fixed. 
Then, we propose to update a according to 



a n+ i = a n - 7? n V a „Var w(a n ). 



(2) 



The sequence {r) n } controls the speed of the adaptive 
process and vanishes as n — > oo, typically like n . The 
novelty of the procedure is that MC sampling and trial 
wave function optimization are coupled. A change in a 
induces a change in the walker dynamical distribution 
which in turn determines the next evolution of a. The 
whole process is non-linear and an explicit numerical in- 
vestigation is required to assess its stability. 

As a specific non-trivial application, we consider the 
two dimensional Gross-Neveu model |l7]] described by the 
Hamiltonian 



H 



dx 



9 



(fV z f) 



(3) 



2N f 

Dirac fermions and we sum over the re- 
N f 



where ip a are Nf 

peated flavor index a = 1, . . . , Nf. The model is asymp- 
totically free, admits a 1/Nf expansion and breaks spon- 
taneously the discrete chiral Z2 symmetry ip — > 75^- 



Following |18| , a lattice formulation with staggered 
Kogut-Susskind fermions jl9| is based on 



H 



L-l 

E 

71 = 



f-(r a ^r a I h r 1 I - <r af r a r at r a Y< 
S 2 ^ n "+ 1 ' gjY^, ^ n C n c n+l c n+l I 



where {c^,c b m } = 0, {c^c^} = S nym 5 a . b and periodic 
boundary conditions are assumed. The state space is 
the set of eigenstates of the occupation number opera- 
tors n° = cf denoted by |n). The fermion number 
is conserved and we focus on the half-filled sector with 
J2i n t = L/2. The Z 2 symmetry corresponds to trans- 
lations by two lattice sites. To avoid sign-problems re- 
lated to boundary crossing we choose in the following 
L mod 4 = 2 (the ground state is then non-degenerate). 
We adopt the one parameter trial wave function 



(n|$(a)) = exp 



L-l N f 



i=0 a=l 



(n|<7 = 0), 



where \g = 0} is the exact ground state at g = 0. 
The algorithm requires an explicit formula for the ra- 
tio (n'|$)/(n|$) where |n) and |n') are states that differ 
by one fermion hopping. If {xi} and {x[} are the L/2 
fermion positions in the two states and if Xi = x[ for 
i =/= p, then the following formula can be derived 



<n'l*) 



e L 



L /2-1 



exp-^ 



exp^p 



Ukjtp ( 



exp-p 



exp — 



We compute the ground state energy on a lattice with 
L = 10 sites and begin our analysis with the case Nf = 2. 
We consider several ensemble sizes and evolution times: 
K = 10, 50, 100 and 500, (3 = 0.1, 0.25, 0.5 and 1.0. 
For each pair (K, f3) we determine by the adaptive algo- 
rithm the best a and estimate the ground state energy. 
For comparison, we also determine Eq by exact Lanczos 
diagonalization. 

Fig. Q shows the typical initial steps of a run. The 
parameter a and the energy measurements evolve and 
fluctuate around (K, 0) dependent definite average val- 
ues a*(K,P) and Eo(K,(3,a*(K,/3)). For large K, the 
statistical error on Eq decreases like K~ x / 2 . For K — » 00, 
the results are expected to be (3 independent. How- 
ever, for moderate ensemble sizes, like those considered 
(K ~ 500), a residual (3 dependence can be observed, 
particularly at intermediate coupling, as shown in Fig. |^. 
This effect is due to the process of walker selection as- 
sociated to SR. The correct approach is to take the 
(3 — > limit where this effect is expected to be negli- 
gible. In Fig. |, we plot E~ (P, 500, a*(500, (3)) as (3 and 
g are varied. All the curves converge to zero and, in fact, 
can be smoothly extrapolated to (3 — > 0. The resulting 
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percentual relative error 100|i?o — Eq\/\Eq\ is very small, 
well below the permille level (see Tab. (I) for numerical 
results with 4th order polynomial extrapolation). 

For large coupling g, the convergence is quite fast. The 
one-parameter trial wave function is accurate because the 
ground state is dominated by states with low potential 
that are easily selected by |$(a)). Relatively small K are 
then already in the asymptotic regime. For intermediate 
couplings, g ~ 2.0, the convergence is again smooth, but 
less than linear. For smaller couplings, a good conver- 
gence is observed and in fact a precise wave function can 
obtained with a* ~ 0. The optimal a* at K = 500, 
(3 = 0.1 is shown in Tab. (I). 

For g = 2.0, we explore a tentative 6-parameter 
trial wave function. Denoting the two fermion fla- 
vors by t, J., we use (n|$) = e^i Fi (n\g — 0} with 

Fi = ELi T~l)+£Lo & *K T 4 +fe + T~l). 

The MC automatic determination of the 6 parameters is 
shown in Fig. ||. The algorithm converges to definite co- 
efficients {a, b}, but the behavior of Eq does not dramat- 
ically improve (Fig. |J). Nonetheless, some qualitative 
remarks can be stressed, as the presence of long range 
correlations between next to neighbor fermions with the 



same spin and anti-correlations between fermions with 
opposite spin. 

Since the Gross-Neveu model can be studied non per- 
turbatively in the framework of the 1 /Nf expansion, it is 
interesting to analyze the algorithm performance with a 
larger number of flavors. In Fig. |(], we show the results 
for Nf = 6. The exact value is beyond Lanczos diagonal- 
ization and we choose to normalize errors at the (3 — 0.1, 
K = 500 value. A comparison with Fig. ||| reveals that 
the error as well as its (3 dependence are rather reduced 
with respect to the previous Nf = 2 case. 

In summary, our data shows that a clever extrapola- 
tion in the algorithm free parameters K and [3 allows 
accurate results even with small walker ensembles. This 
is an important feature for realistic large scale simula- 
tions aimed at reaching the continuum limit. Results 
with large Nf suggest that the present algorithm can 
be a viable numerical technique for other fermionic two- 
dimensional models where the 1/Nf expansion applies, 
like the important case of models with dynamical super- 
symmetric breaking |2(J. In principle, extensions to mod- 
els with sign-problems are possible and, in fact, progress 
in the optimization issue has been recently proposed pl| ] 
within the considered class of MC algorithms. 



TABLE I. Eo/N f for the L = 10 model with N f = 2 flavors. AE = E% anczos - E ( 



g 


a* (500, 0.1) 


Exact Lanczos Diagonalization 


Polynomial Extrapolation 


1000 \AE/E\ 


0.5 


0.07638(1) 


-3.34904 


-3.34908(5) 


0.012 


1.0 


0.31347(5) 


-3.71687 


-3.71689(5) 


0.005 


2.0 


1.4044(3) 


-5.99265 


-5.9929(5) 


0.03 


2.5 


2.0575(2) 


-8.4526 


-8.4524(3) 


0.02 


3.0 


2.6198(2) 


-11.6949 


-11.6927(3) 
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FIG. 3. L = 10, N f = 2, K = 500. Relative percentual 
error on the energy obtained from data at large K at several 
/3. 
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Monte Carlo iteration 



Monte Carlo iteration 

FIG. 1. L = 10, N f = 2, g = 3.0, K = 10, /3 = 0.5. MC 
evolution of the ground state energy estimate and of the a 
parameter. 



FIG. 4. L = 10, N f = 2, g = 2.0. MC evolution of the six 
parameters {a, b}. From top to bottom, on the right of the 
plot, the parameters are a\, a^, ai, bi, b%, bo- 



g = 0.5 



< 1 






g = 2.5 
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FIG. 2. L — 10, Nf = 2. Relative percentual error on the 
ground state energy. The various lines correspond to j3 = 0.1 
(circles), /3 = 0.25 (squares), f3 = 0.5 (diamonds), /3 = 0.75 
(triangles up) and /3 = 1.0 (triangles down). 
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FIG. 5. L = 10, Nf = 2, g = 2.0, K = 500. Improve- 
ment in the energy estimate with the 6-parameter trial wave 
function. 



< 0.5 
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FIG. 6. L = 10, N f = 6, K = 500. Relative percentual 
error on the energy estimate obtained from data at large K 
at several f3. 
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